Bose-Einstein condensates with a bent vortex in rotating traps 
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We consider a 3D dilute Bose-Einstein condensate at thermal equilibrium in a rotating harmonic 
trap. The condensate wavefunction is a local minimum of the Gross-Pitaevskii energy functional 
and we determine it numerically with the very efficient conjugate gradient method. For single vortex 
configurations in a cigar-shaped harmonic trap we find that the vortex line is bent, in agreement 
with the numerical prediction of Garcia- RipoU and Perez-Garcia, Phys. Rev. A 63, 041603 (2001). 
We derive a simple energy functional for the vortex line in a cigar-shaped condensate which allows 
to understand physically why the vortex line bends and to predict analytically the minimal rotation 
frequency required to stabilize the bent vortex line. This analytical prediction is in excellent agree- 
ment with the numerical results. It also allows to find in a simple way a saddle point of the energy, 
where the vortex line is in a stationary configuration in the rotating frame but not a local minimum 
of energy. Finally we investigate numerically the effect of thermal fluctuations on the vortex line 
for a condensate with a straight vortex: we can predict what happens in a single realization of the 
experiment by a Monte Carlo sampling of an atomic field quasi-distribution function of the density 
operator of the gas at thermal equilibrium in the Bogoliubov approximation. 

PACS numbers: 03.75.Fi, 67.40.Vs 



I. INTRODUCTION 

Several experimental groups have now produced vor- 
tices in Bose condensates of atomic gases by a rotation of 
the trapping potential ^. It is therefore important 
to characterize the steady state vortex configurations for 
parameters relevant to the experiments. In particular 
the regime of the ENS and MIT still deserves some in- 
vestigation. It corresponds to a condensate trapped in a 
cigar-shaped harmonic potential, that is with a harmonic 
oscillation frequency along the rotation axis z smaller by 
one order of magnitude than the (quasi-degenerate) os- 
cillation frequencies in the x — y plane. Furthermore, a 
configuration with a single vortex can be produced in the 
experiment in a reproducible way The single vortex 
configuration in a cigar-shaped trap is the main object 
of the present work. 

Why is the single vortex cigar-shaped regime so in- 
teresting? First the weak harmonic confinement along z 
makes the system very different from previously studied 
configurations. If the harmonic confinement along z was 
stronger than in the x — y plane the vortex configura- 
tions would be similar to what happens in a 2D rotating 
Bose e well studied limiting case [|[ ^, ||, ^ . If the 
harmonic confinement along z was absent, e.g. if the con- 
densate was a cylinder, one would face a situation typical 
of superfluid hehum II, a well studied subject [^j. 

But now the gas keeps its 3D character while being 



spatially inhomogeneous along z. This has the amusing 
consequence that a single vortex line has a tendency to 
bend. First the steady state condensate can have a bent 
vortex line, as shown numerically by a minimization of 
the Gross-Pitaevskii energy functional in Q and as ob- 
tained also numerically with an approximate vortex line 
energy functional in pO[ . Second, even a condensate with 
a straight vortex line has low energy Bogoliubov modes 
corresponding to large fluctuations of the end points of 
the vortex line [^l], |l2| so that the vortex line can easily 
bend in presence of thermal fluctuations. 

After a presentation of the model considered in this ar- 
ticle, see section ^ we give a systematic numerical study 
of the steady state bending as function of the trap as- 
pect ratio, see section ID, an analytic understanding of 
the bending of the steady state vortex line in a cigar- 
shaped trap, see section IV, and the description of the 



finite temperature fluctuations of the vortex line, see sec- 
tion including a discussion of the effect of vortex line 
bending on the experimental absorption images. 



II. MODEL AND BASIC ASSUMPTIONS 

Let us consider an almost pure Bose-Einstein conden- 
sate of N atoms confined in a trapping potential rotating 
along axis z at angular velocity Q. The thcrmodynam- 
ically metastable equilibrium configurations of the sys- 
tem correspond to local minima of the Gross-Pitaevskii 
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energy functional in the rotating frame 
E\<h,ch*] = / £r 



{Ho-nL.)<l>+^\<l>f 



(1) 



where the condensate wave function (f) obeys the normal- 
ization condition 
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(2) 



The Hamiltonian operator Hq in Eq. contains the 
kinetic and trapping potential terms 
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U{r) 



(3) 



and Lz — ^ifi{xdy — ydx) is the angular momentum op- 
erator along z. Here we will consider harmonic trapping 
potentials U{r) with an adjustable slight anisotropic de- 
formation in the x ~ y plane: 

U{f) = ^mul [(1 - e)x^ + (1 + e)y^] + ^mco^z^. (4) 

For the choice of parameters we refer to the typical values 
of the recent ENS experiments with ^^Rb atoms (scatter- 
ing length a = 100 oq — 5.29 x lO^^m) [Q: total number 
of atoms N = 1.4 x 10^, axial frequency = 27r x 11.7 Hz 
and anisotropy parameter e 1. In this paper the radial 
frequency uj± is varied in order to explore a wide range 
of trap anisotropics: = uj±/iVz € [1, 15]. In what fol- 
lows en ergies an d lengths are given in units of huj± and 
a_L = y^h/mLu± respectively, m being the atomic mass of 
87Rb. 



III. LOCAL MINIMA OF ENERGY: 
NUMERICAL RESULTS 

In this section we discuss the method and the algo- 
rithm employed to minimize numerically the energy func- 
tional (|^). Then we present the results obtained for the 
stationary configurations with and without vortices in a 
long cigar trap, and finally we discuss the existence do- 
main for a single vortex configuration for several trap 
geometries. 



A. Method 

First of all, we reformulate the energy functional in 
order to account more easily for the normalization con- 
straint (H): we define 4> = 4}/\\tp\\ so that the value of 
E can be obtained for condensate wave functions ip with 
a norm different from unity. This corresponds to divid- 
ing the terms of E[tp,^*] quadratic in tp by and 
the interaction term quartic in tjj by HV"!!^- The modified 
energy functional reads 



Ng IVI^ 



Then we discretize ^ on a, three-dimensional grid with 
periodic boundary conditions in position and in momen- 
tum space. The number of grid points that we use ranges 
from 64 to 256 for each lattice side, depending on the 
trap geometry (a typical choice is 64 x 64 x 192 for long 
cigar configurations and 96 x 96 x 64 for spherical ge- 
ometries with vortices). The minimization is performed 
by using the conjugate gradient algorithm described in 
p3[ . This algorithm is in general much more efficient 
that the usual steepest descent method (see Ref. jl^ 
for a comprehensive discussion). One ingredient of the 
conjugate gradient method is a line minimization of the 
energy functional, that is the minimization of E along 
the line ^ — ipo + XSip, A being a real parameter, where 
tAo is the current trial wave function and S'tjj is the pro- 
posed direction (or gradient) along which to move the 
trial wave function. An important issue here is to find 
the first minimum encountered when moving downhill in 
energy along the line: if the algorithm can jump to an- 
other minimum on the line, corresponding for example to 
an energy valley with a different number of vortices, one 
gets wrong predictions, in the sense that one does not 
fully explore the stability domain of a branch of solution 
with a given number of vortices. This issue is usually not 
considered as important in the textbook implementation 
of the conjugate gradient method: e.g., in ||l3| the line 
minimization does not necessarily find the first minimum. 
What we have done is to use the fact that £^ is a rational 
function of A so that it is completely characterized by the 
coefficients of the monomials in A in the numerator and 
the denominator. We then easily find the roots of dE/dX 
and the first local minimum of E encountered when one 
moves along the line downhill in E starting from A = 0. 
Afterwards we proceed with another line minimization 
along a conjugate direction p^ , and so forth until we 
find a local minimum of the energy functional (|^) in the 
full configuration space spanned by the wave function "0. 
As convergence criterion we use 



\\Hgpip - ntPW < 



where Hgp is the Gross-Pitaevskii Hamiltonian 



and fi is given by 
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(5) 



(6) 
(7) 

(8) 



which eventually converges to the chemical potential. 
The parameter 7 is a small parameter (~ 10^® 10^^") 
controlling the quality of the convergence. 



B. Results for a fixed trap geometry: case of a long 

cigar 

We start by considering a cigar shaped trap geometry 
typical of the ENS experiments in which vortex configu- 
rations are nucleated [0, |ll|. In particular we set e = 0.03 
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and uj± = 2tt X 174Hz, giving ~ 14.9. To explore 
the configuration space we start from a trial wave func- 
tion for a fixed value of the rotation frequency fl until 
an equilibrium configuration (i.e. a local minimum of 
the energy functional (|^)) is found. Then we increase 
or decrease slightly f2, and then find the new local min- 
imum. In this way we can follow continuously branches 
of configurations with or without vortices, depending on 
the rotation frequency and on the path followed. In Fig- 
ure |l| we show the energy E and the angular momentum 
per particle (L^) along the rotation axis of each config- 
uration, as function of the rotation frequency Q, up to 
4-vortex configurations. 

We notice that, for this trap geometry, when a single 
vortex appears (i.e. when it becomes a stable config- 
uration) it has immediately a lower energy than the 0- 
vortex configuration, differently to what happens in 2D 
1^. Moreover this vortex configuration is characterized 
by a bending of the core line, a phenomenon already ob- 
tained with a different numerical technique in and 
this accounts for the fact that the angular momentum 
per particle (Lz) is lower than h {{Lz) is equal to fi for a 
straight centered vortex). To better investigate this as- 
pect we have studied the deformation of the vortex line 
as function of the rotating frequency. If we start with a 
single vortex configuration and we decrease the rotation 
frequency 12 the bending of the core line becomes more 
pronounced, until we reach a certain critical frequency 
r^i and the system jumps to the 0-vortex configuration. 
In the opposite direction, when il increases, we find that 
the vortex line tends to become more straight as expected 
(see Figure^), but then, whereas a small bending is still 
present at the extremities of the condensate, it "decays" 
into a three vortex configuration. The value of O at which 
this happens defines a second critical rotation frequency 

Therefore for the long cigar trap and the parameters 
considered here, one never obtains a single, straight vor- 
tex: when f2 is too large, other vortices come in. As we 
will show in the following, this is due to destabilization 
of high angular momentum surface modes [p^ , and cor- 
responds to the fact that the single vortex configuration 
ceases to be a local minimum of energy when the energy 
of these modes becomes negative (i.e. there is at least 
one direction along which the energy is going down in the 
functional space). We discuss this effect in details in the 
next subsection. 



C. Destabilization of surface modes 

We consider in this subsection the case of a cylin- 
drically symmetric trapping potential (with a vanishing 
asymmetry parameter e = 0) so that eigenstates of the 
A^-body Hamiltonian have a well defined angular momen- 
tum. The expression for the energy £i of an excitation 
mode of the condensate with angular momentum Til in 
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FIG. 1; Energy per particle (top) and angular momentum per 
particle (Lz) (bottom) as a function of the rotation frequency 
57 for various branches. Circles: starting from the 51 = 
ground state without vortices we reach configurations with 
2 and 4 vortices as we increase 57. Triangles: configurations 
with one bent vortex which decays into a 3 vortex (no vortex) 
state as we increase (decrease) 51. This branch is obtained by 
starting from a trial wave function with quantized circulation 
of charge -1-1 for a value of 57 which can support a vortex state. 
See the text for the value of the parameters. The energy is in 
units of Tiijj± and (Lz) in units of Ti. 



the rotating frame is 

£i = huj{i) - inn, 



(9) 



where uj{l) is the excitation energy of the same mode in 
the laboratory frame This determines the rotation 
frequency fi; = toil) /I at which the energy of such mode 
becomes negative in the rotating frame. In particular the 
minimum value over I of fi; defines the Landau critical 
rotation frequency f2 at which the vortex becomes ther- 
modynamically unstable due to destabilization of these 
surface modes 
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(10) 



The value can be estimated analytically by using the 
sum rules approach which provides an estimate for uj(l), 
as discussed in Refs. ||lj, |l^. In order to describe the 
surface modes, we consider the excitation operators for 
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Then one defines the moments 
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±^ / dE[S+{E)±S^iE)]EP 
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(12) 



where the strength distribution functions S± are given 
by 



S±iE) = J2\{n\F±\0)\^SiE - hcj,,). 



(13) 



The states |n) form a complete set of eigenstates of the 
Hamihonian operator H for our system of N interacting 
particles confined by the trapping potential: 



N 



fc=i 
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fe-i 



(14) 



Here we assume the two mode approximation, as dis- 
cussed in Rcf. 
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FIG. 2: Density cuts in y — z plane (vertical and horizontal 
axes respectively) of a long cigar condensate (A"'^ — 14.9) 
with a bent vortex, for three values of the rotation frequencies: 
Q/lu± — 0.37,0.44,0.59 from top to bottom. Lengths are 
given in unit of a±. 
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where the two modes are equally weighted (they have the 
same strength a) due to the vanishing of the rriQ momen- 
tum This is justified by the fact that in the large 
N limit the strength distributions are dominated by the 
contribution of two modes with energy 'hLij± and angu- 
lar momentum ±?i/ excited respectively by the operators 
F±. With this ansatz it is straightforward to obtain the 
relation between the frequencies uj± of such modes and 
the first three moments in Eq. (|l^ ) ; in particular we have 
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UJj^ — LO^ = 
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(16) 
(17) 



Notice that only the mode with frequency w+ = a) -I- 5 is 
relevant to our discussion, since the energy of the mode 
with angular momentum —hi increases for increasing ro- 
tation frequency f2, and can never become negative. 

The next step is to calculate the moments from 
Eqs. (|lj) and (|l^) by using closure relations and then 
explicitly evaluating the commutators 
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([F_,[i7,^^+]]) = ^/2(rl'-) 



= ([[F_,i7],[i?,F+]]) = 
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(18) 
(19) 
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2h{l - 2){rf-^L,) + 2(1 - 2)(rl'-^L2)/? 
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where (. . .) stands for the expectatfon value in the state |0). These expressions generahze those in (case of a 
vortex free condensate) and in [|l5) (case of excitations of Z = 2 of a condensate with vortices) . We note that the use 
of p± = px^ipy and r± = x±iy rather than the Cartesian coordinates Pa: , Py, x, y, e.g. by writing the kinetic energy 
as p+p-^ considerably simplifies the calculations of the commutators. 

For the case of a straight vortex with charge q the key quantities to be inserted in Eqs. ( p^ and ( p7| ) are 
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P-) , 2q{l + q){l^2){l-l){rl-^) 
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(21) 
(22) 



where quantities within square brackets are expressed in 
dimensionless units, as defined in Sect. ||. 

In Figure ^ we plot the curve ui[l)/l for a cylindri- 
cally symmetric trap (e = 0), for the case q = Q (ground 
state without vortices) and g = +1 (straight vortex). 
In the same figure we also indicate the critical frequen- 
cies ^2 ^2 which define the upper bound of the ex- 
istence domain respectively for states without vortices 
and with a bent vortex, as found from the numerical 
minimization of the energy functional in presence of a 
slight anisotropy as considered in the previous subsec- 
tion (e — 0.03). We notice that the sum rule prediction 
gives a value of Toim{uj{l) /I) which is quite close to our 
numerical result. This supports the statement that the 
surface modes are indeed responsible for the destabiliza- 
tion of the single-vortex state and the vortex-free state. 
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D. Effect of a ballistic expansion 

In general, in the experiments, the condensate density 
is imaged after a time-of- flight: the trapping potential is 
switched off and the cloud expands freely for some time 
before being imaged by the absorption of a laser beam. 
This has the advantage of increasing the diameter of the 
vortex core, making it larger than the optical wavelength. 
Some theoretical work is however needed to relate the 
images obtained after the expansion to the density profile 
of the initial trapped condensate. 

To calculate the expansion of the condensate after the 
release from the trap, we propagate the initial wavefunc- 
tion according to the full 3D Gross-Pitaevskii equation 
in real time after having abruptly switched off the trap- 
ping potential, and after having performed the scaling 
and gauge transform of . In terms of the rescaled 
spatial coordinates fi = ri/Xi{t), the modified wave func- 



FIG. 3: Curve lu{1)/1 for a condensate in a long cigar trap 
(A~^ = 14.4, e = 0) with (squares) and without (circles) a 
strai ght vortex, as calculated with the sum rules approach (see 
Eq. ( |2l| ) and (|22[)). The minimum of uj{l)/l, which gives an 
estimate for the rotation frequency Q2, is compared with its 
corresponding value found numerically for the 0-vortex (f22i 
dot-dashed line) and the 1-vortex {^2, dashed line) states, in 
presence of a slight anisotropy (e = 0.03). Frequencies are 
given in units of uj±. 



tion ip{r, t) satisfies the equation 
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and the scaling parameters \j (t) are solutions of ||Tf 



A, 



A> A1A2A3 



(24) 



with initial values equal to unity and with vanishing ini- 
tial derivatives. The use of these scaling equations allows 
us to calculate numerically the expansion of any vortex 
configuration by using the same grid as the one used for 
the stationary trapped state as most of the evolution of 
■0 due to the ballistic expansion is absorbed by the scal- 
ing and gauge transform. In a 2D model with an axially 
symmetric trapping potential the net effect of the bal- 
listic expansion on the density corresponds exactly to a 
scaling transformation and a finite angle rotation of the 
condensate eigenaxes Here we expect this property 
to remain approximately true as the condensate is very 
elongated along z and expands much more rapidly radi- 
ally than axially. 

A first application of this calculation is the analysis of 
the column density of the condensate, that is the den- 
sity integrated along z, after ballistic expansion. There 
is indeed an important issue concerning the mechanisms 
which are responsible for the reduction of contrast for 
the vortex configurations which are observed in the ex- 
periments As pointed out already in in the case 
of a bent vortex before ballistic expansion, an important 
contribution can arise due to the bending itself, as shown 
in Figure ^ In this picture we plot a cut of the column 
density along the y axis where the vortex bends, for three 
values of the rotation frequency f2. In Figure g we com- 
pare the column density for a bent vortex ai fl ~ 0.59u!±, 
before and after the release from the trap (time of fiight 
~ 30ms). The y coordinate is given in rescaled units 
y = y/\2{t). From this figure we can see that the ex- 
pansion has no relevant effects on the contrast; the same 
result holds also for the other configurations shown in 
Figure ||. 

A second application is the analysis of the shape of the 
vortex line after ballistic expansion. This is also an im- 
portant issue, as the shape of the vortex line can be mea- 
sured using transverse absorption imaging rather than 
axial absorption imaging p^ . Focusing on the configu- 
ration O = 0.37a;_L exhibiting the strongest bending, we 
have calculated the shape of the vortex line before and 
after expansion by looking in each horizontal plane for 
the local minimum of density closest to the vertical axis 
z. Initially the vortex line is contained in the half plane 
{—y)z which is at an angle 6{Q) = — 7r/2 with respect 
to X axis. After a time — 30ms of ballistic expan- 
sion, we find that the vortex line remains almost planar, 
but in a vertical half plane at an angle 9{te) — — 1.14rad 
with respect to x axis, and not passing through the ori- 
gin. We show in Figure |^ the vortex line before and after 
expansion, projected in the vertical planes in which it is 
contained at f = and t = t^- One then sees that the vor- 
tex line essentially preserves its initial shape after time 
of fiight, in the frame of the rescaled spatial coordinates. 
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FIG. 4: Column density before ballistic expansion (inte- 
grated density along the axial direction) along the y axis 
where the vortex bends, for the three cases in Figure ^: 
Q,/uj± = 0.37, 0.44, 0.59. This picture evidences that the con- 
tribution of the vortex line bending to the contrast can be 
relevant even at zero temperature. 
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FIG. 5: Comparison of the column density for a bent vortex 
{Q, = 0.59(j;x) before and after the release from the trap (time 
of flight ~ 30ms), as a function of the rescaled y coordinate. 
The expansion has no relevant effects on the contrast (we have 
verified that this holds also for the other configurations shown 



in Figure ^) 



E. Existence domain of a single vortex 
configuration 

In the last part of this section we discuss the thermo- 
dynamic existence domain of a single vortex configura- 
tion, that is the domain where it is a local minimum of 
energy. By varying the radial frequency w_l we have ex- 
plored a wide region of trap geometries, ranging from the 
spherical case to long cigar traps (A^^ G [1, 15]). In Fig- 
ure 1^ we show this domain in the u!±/lOz — ^/u>± plane. 
Here we consider an almost axially symmetric trap, with 
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FIG. 6: Shape of the vortex hne (with initial rotation fre- 
quency Q = 0.37tJ_L) before and after release from the trap 
(time of flight te ~ 30ms), as a function of the rescaled co- 
ordinates in units of a±_ and after projection on the vertical 
planes at an angle 6(0) = — 7r/2 and 6{te) — — 1.14rad with x 
axis, respectively. 



a very small asymmetry (e = 10~*), which fixes the di- 
rection along which the vortex bends. By minimizing 
numerically the energy functional (|^) with the conjugate 
gradient method we find that thermodynamically stable 
configurations with one vortex lie in the interval between 
the lines with empty and full circles, representing respec- 
tively f2i and Q2- The black squares represent the pre- 
dictions of the sum rule s, wh ich give a good estimate 
for ri2, as discussed in jjlllC, 



In the same picture we 
also show the results from the Bogoliubov approach for 
a straight vortex, as described in §0. Up-triangles cor- 
respond to the critical rotation frequency which stabilize 
a straight vortex (f2i); this line separates the existence 
domains for straight and bent vortices. Down-triangles 
correspond to the frequency at which the straight vortex 
is no longer thermodynamically stable due to the desta- 
bilization of the surface modes (O2). Notice that these 
values are in very good agreement with the result of the 
conjugate gradient in the whole range of existence of the 
straight vortex. 

Finally, the solid line represent a 2D ansatz based on 
the analytical model in Ref . |^ , as described in the next 
section. If we imagine our 3D condensate as a collection 
of 2D slices, we can suppose that the 3D vortex is a sta- 
ble configuration {i.e. the vortex line remains close to z 
axis and does not get away) if the rotation frequency 
at least exceeds the 2D critical rotation frequency 
at which a 2D vortex in the central slice {z — 0) becomes 
energetically favorable with respect to the solution with- 
out vortices, that is: SI > ^l^{z = 0). The expression 
for fi^^, which will be derived in the next section, is 
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FIG. 7: Phase diagram for the existence domain of a single 
vortex as a local minimum of energy in an almost cylindrically 
symmetric trap. For a given A, thermodynamically stable con- 
figurations lie in the interval [r2i,r22]. The Bogoliubov pre- 
dictions from fli and il2 (triangles) correspond to a straight 
vortex. In the conjugate gradient minimization (disks) the 
bending of the vortex line is allowed and the stability domain 
is enhanced. 



where the tilda indicates rescaling by e.g. = 
^i/hLj_i_, C ~ 1.8011, and rj = 1.9378. Then, by ex- 
pressing /i in terms of A = lUz/lo^ in the Thomas-Fermi 
approximation: 



we have 



15iV- 



2/5 
Ai/5 



13. lA^/^ 



Ai/5 



In ( cA^/-"^ + d 



(26) 



(27) 



with b ~ 0.0763, c ~ 23.6 and d ~ 3.49. As we can see 
from Figure this simple ansatz gives a good estimate of 
the critical frequency Qi. This ansatz will be justified in 
the next section. 



IV. UNDERSTANDING THE BENDING OF 
THE VORTEX LINE ANALYTICALLY 

The previous numerical results show that the vortex 
line in a steady state configuration is not necessarily 
straight when the condensate is cigar shaped along the 
rotation axis z, in accordance with previous numerical 
results based on a different algorithm This how- 
ever does not explain physically why the vortex bends. 
To get the required physical understanding we derive an 
approximate energy functional for the vortex line in the 
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Thomas- Fermi limit in the spirit of |L(][ , and we minimize 
this energy functional with a simple variational ansatz. 
We reach a very simple prediction for the minimal rota- 
tion frequency required to stabilize a bent vortex, which 
is in good agreement with the full numerics when the 
condensate is cigar-shaped. 



A. Deriving a simple energy functional 

We restrict in what follows to the interesting regime 
of a cigar shaped condensate, where lOz <C wj^. The first 
step is to transform the Gross-Pitaevskii energy func- 
tional (0) into a functional of the shape of the vortex line 
only. This assumes that both the modulus and the phase 
of the condensate wavefunction can be expressed in 
terms of the shape of the vortex line only. This is in gen- 
eral a formidable task, as the condensate density is not 
uniform in a harmonic potential ]l9| , but it is greatly sim- 
plified if we restrict to the Thomas- Fermi limit /i ^ fiLij±^ . 
We present a rather detailed and pedestrian derivation in 
the appendix |a|, we give here only the main ideas. 

In the Thomas-Fermi regime there is a clear separa- 
tion of spatial scales between the vortex core radius, of 
the order of the healing length ^ — (?i^/m/j,)^/^, and 
the transverse Thomas-Fermi radius of the condensate 
R±_ = i2ijL/mu)Y)^^'^ , where n is the chemical potential 
of the gas. The total density can then be written as the 
product of a slowly varying envelope function and of a 
narrow 'hole' function defining the vortex core |6|. We 
further assume that the rotation frequency il is of the 
order of hu>\/ ji. As a consequence the rotational veloc- 
ity term 17 A r, at most of the order of r2i?_L, is much 
smaller than the typical velocity field in the lab frame 
at a distance ^ from the vortex core, v ~ h/m£,, in the 
Thomas-Fermi limit, so that the structure of the vortex 
line is not distorted by rotation. Another consequence 
is that the envelope function is also not destabilized by 
the rotation and is close to the usual Thomas-Fermi 
expression. 

Expressing the phase of the condensate as function of 
the vortex line is made difficult by the spatial inhomo- 
geneity of the density profile of the condensate In 
principle this phase has to be determined everywhere in 
the condensate if one wants to calculate the kinetic en- 
ergy term of (|^) . Fortunately, using the fact that the con- 
densate is in a steady state, one can replace the volume 
integral giving the kinetic energy stored in the conden- 
sate phase by a line integral along the vortex line, using 
the same type of techniques as in 0. It is then possi- 
ble to rely on approximations for the condensate phase 
valid close to the vortex line. We use in the appendix ^ 
the simplifying hypothesis that the vortex line is weakly 
curved, with a radius of curvature of the order of R± or 
larger, which allows to approximate the condensate phase 
close to the vortex line by the one of a straight vortex. 

We finally obtain the following energy functional of 
the vortex line, taking the vortex free configuration as 



the zero of energy: 



ds 
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+ cos{a{s))E:^'^{ro±{s);z„{s))] . (28) 



In this expression the vortex line is parametrized by the 
curvilinear abscissa s. At the point of abscissa s the 
vortex line is at the elevation Zq{s) and at a distance 
^O-l('S) from the rotation axis, and makes an angle a{s) 
with respect to z. A remarkable feature of (|2^) is that 
it is expressed in terms of the energy functionals of a 
vortex core in a 2D condensate, for the energy in 

the absence of rotation and i?2D ^oi the energy due to the 
—^ILz term. This is physically plausible considering the 
cigar shaped nature of the condensate, and this allows to 
view the 3D condensate as a collection of 2D horizontal 
slices. The slice of elevation z constitutes a 2D Bose 



condensate with a chemical potential fj,- 



^2, where 



/i is the chemical potential of the 3D condensate, and 
has a Thomas-Fermi radius coinciding with the local 3D 
one. The coupling constant (/2d(-z) of the 2D gas can 
be ex pressed in terms of the 3D coupling constant, see 
( A58 ) . We arrive at the simple formula 
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15 



2 \ ' 



(29) 



where Rz is the Thomas-Fermi radius of the condensate 
along z. An interesting remark is that the rotational 
energy term in ( |2^ ) is multiplied by cosa(s). As £'213 
is proportional to the rotation frequency Q, this means 
that cos q;(s)£'2d is the rotational energy of a vortex core 
in a 2D condensate rotating at the effective frequency 



f^2D(s) = ricosa(s). 



(30) 



Minimizing numerically the simple energy 
functional 



In a first stage we have to check that the energy func- 
tional derived in the appendix ^ correctly reproduces the 
results of the minimization of the full Gross-Pitaevskii 
energy functional. We perform this check numerically: 
we discretize the vortex line in little segments having 
all the same length dl much smaller than the transverse 
Thomas-Fermi radius R± of the condensate. As the vor- 
tex line is symmetric with respect to z reflexion and lies 
in the xz plane, the left extremity of the first segment in 
the calculation moves along x axis only, with an abscissa 
xq. The z > part of the vortex line is discretized in k 
segments and its shape is parametrized by the angles a^, 
I = 1, . . . , fc at which the k segments are with respect to 
the axis z. The energy functional (|2|) in its discretized 
version is then a function of -I- 1 coordinates, that is of 
Xq and of the k angles ai . Starting with a straight vortex 
line at some small angle with respect to z axis we move 
Xq and the ai 's according to the simple gradient method 
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FIG. 8: Minimal rotation frequency required to stabilize a 
vortex in a cigar-shaped condensate, as function of A^^ = 
tJx/<^z- The parameters are given in the last paragraph of 
section ^ Disks with error bars: minimization of the full 
Gross-Pitaevskii energy functional. Squares: minimization of 
the approximate vortex line ener gv fu nctional, based on the 
2D 77- modified energy functional (4.66) and a discretization of 
the z > part of the vortex line in k — 256 segments. Solid 



line: analytical estimate of subsection IV C 



or imaginary time evolution method, that is we move the 
parameters by a small step in the direction opposite to 
the local gradient of the energy functional. 

It is known that this simple gradient method is not 
efficient when the desired minima are at the bottom of a 
very elongated valley in the coordinate space This 
potential problem is minimized by a rescaling of the coor- 



dinates by their natural units x, 



typ 



i? I and a*yP = 2tt, 



so that our specific algorithm is to iterate the following 
small coordinate changes 



dxo - ~dT{x'yP) d^„E^/E'JP (31) 
da, = -dT{a'ypYda^E,/ElyP (32) 



where the typical energy scale is E^^ — fi^uj^/ ^ and the 
dimensionless 'time' step dr is small as compared to unity 
(typically 0.1). The iteration stops when \dxQ\/ {drx^^^) 
and \dai\/ {dra^^^) are below some threshold, taken here 
to be 10"^. 

The resulting prediction for the minimum rotation fre- 
quency required to stabilize the vortex line is shown in 
Figure ||, deliberately restricted to the domain of cigar 
shaped condensates. The agreement with the minimiza- 
tion of the full Gross-Pitaevskii energy functional is re- 
markable, considering the fact that the points of Figure^ 
are moderately in the Thomas-Fermi regime (/i ~ Ldhajj^ 
for uj±/ujz = 15). 



Why the vortex line bends in a cigar shaped 
condensate 



The actual goal of this section is to understand phys- 
ically why the vortex line bends in a cigar shaped con- 
densate. This can be achieved intuitively thanks to the 
very suggestive form of the energy functional (28). One 



just needs to have in mind the following characteristics 
of the 2D vortex problem: 

• if the effective rotation frequency f22D is too small 
the 2D energy functional has a maximum for the 
vortex core at the center of the trap and is a purely 
repulsive potential, see the dotted line in figure ||: 
the vortex core cannot be stabilized inside the con- 
densate and its equilibrium position is at infinity. 

• if ri2D is above the stabilization frequency ^^^^^ 
and below the critical rotation frequency ^c^, the 
2D energy functional has a local but not global min- 
imum for the vortex core at the center of the trap 
(dashed line in Figure ||). In this situation, the 
vortex core is stabilized at the trap center with an 
energy larger than the vortex free condensate. 

• for f22D > ^^c'^, the energy minimum at the trap 
center is now below the energy of the vortex free 
configuration (solid line in Figure |^) . 

The important feature of the 2D case is that the equi- 
librium positions of the vortex core are either the trap 
center or the infinity. Another point, crucial for the 3D 
case, is that both rigt'ab ''^^'^ ^""^ decreasing functions 
of the chemical potential /Lt2D- 

Let us now follow the vortex line travelling through 
the 3D cigar shaped condensate, in the case where the 
2D critical rotation frequency in the central slice z — 
is smaller than the actual rotation frequency f2. Let us 
call Zc the elevation of the 2D slice with a local critical 
frequency fi^^ equal to see Figure [lC|a. 

It is clear that the vortex line will be straight along 
the rotation axis for \z\ < Zc'. the vortex line is there 
in a valley corresponding to the global minimum of en- 
ergy of each local 2D slice. When the vortex line reaches 
the domain of elevation z > Zc, having a vortex core on 
the rotation axis costs more energy than having the vor- 
tex core at infinity in each local 2D slice. The tempting 
strategy then offered to the vortex line is to bend and 
leave the condensate. Assume that the vortex line leaves 
the condensate radially, as shown in Figure The 
corresponding horizontal vortex line has an energy 



E^ 



horiz 



92-d(zc 
9 



+ 00 



dxE^^"{x;z 



(33) 



The corresponding integral can be calculated exactly, but 
it is sufficient to give here an order of magnitude: 920/9 
scales as l/i?z, -E^ff" is of the order of h^u)j_/ ^21) and 
the integral over x converges over a distance given by the 
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FIG. 9: Energy of a single vortex in a 2D Thomas- Fermi 
condensate in a harmonic trap mLu'^{x'^ + J/^)/2 rotating at 
frequency fi, as function of the distance of the vortex core 
to the trap center, for a chemical potential fj,2o ~ 10hu!±. 
Dotted line: Q — 0.1lj±. Dashed line: Q — 0.25lj±. Solid 
line: Q = 0.35tJx. The values of the stabiliz ation and critical 
rotation frequencies defined in the text, see ( |A7C| ) and (A68), 
are Q^°^^ ~ 0.170lo± and 0,'^° ~ 0.306cj±. The energy is 
in units of h^cu'^/ (^20 where /i2D is the chemical potential of 
t he 21) gas, and is calculated from the approximate formulas 
( A64 , A65 A66 ) . The distance is in units of the Thomas- Fermi 
radius of the 2D condensate. 



Thomas- Fermi radius so that 



E. 



horiz 



3/2 



(34) 



We have included some approximate Zc dependence rele- 
vant for the extreme case of Zc close to Rz ■ What would 
be the energy cost for the vortex line to remain on the 
rotation axis from Zc to i?^? The energy of the corre- 
sponding vertical segment is 




FIG. 10: The two C-shaped ansatz used in this paper for the 
vortex line in a cigar-shaped condensate in a trap rotating at 
frequency Q,. (a) the vortex line (thick solid line) is on the 
rotation axis for an elevation in between —Zc and Zc, otherwise 
goes to infinity horizontally; Zc is the positive elevation of the 
2D slice having a critical rotation frequency £7^° equal to f2. 
(b) the C shape is at a distance X to the rotation axis and 
has a total height 2Z\ X and Z are determined variationally. 



the strategy of bending is more favorable than the strat- 
egy of following the rotation axis, except for a Zc very 
close to the end point of the condensate, 1 — z^jRz ~ 
R\lRz that is r2 ~ {r?up-J pi){Rzl Ri^"^ . 

The bending at the value z^ in the above reasoning 
can be justified variationally as follows. We perform a 
simple minimization of the energy functional ( psf ) with 
the following linear piecewise variational ansatz for the 
vortex line in the plane y = 0: an horizontal line linking 
(x — -|-cx), z = — ^) to (x = 0, z = — ^), then a vertical 
segment of length 2Z along the rotation axis, then an 
horizontal line from (x = 0, 2: = Z) to (x = -|-c», z = Z). 
The energy of the ansatz depends on the single varia- 
tional parameter Z: 



EZ 



dz ■ 



15 
16E 



[E^^\0-z)+El^{0-z)\ (35) 

""'dz {i-^^\nnl^{z)-n^l''{z^)\ 



where we have used the fact that the 2D rotational energy 
of a vortex core in the center of a trap rotating at fre- 
quency ri is —Ml, also equal to —Ml1^[zc) by definition 
of Zc- The integrand has the same order of magnitude as 
in E^°™ but the integration length is now of the order 
of Rz so that 



EZ 



Zr 
1 

Rz 



(36) 



is typically Rz/R± times larger than E^""^^. This proves 
that in the limit of a cigar-shaped condensate Rz S> R± , 



E^'''{Z) = sm^^') + E^^'^'^iZ) 



(37) 



that is the sum of the energies of the vertical segment 
and of the horizontal lines. One has to extremize this 
function over Z. In the case where the extremum is in 
the interior of the interval (0, Rz), one has to solve 



dZ 



ET^ = 0. 



(38) 



This non-trivial task becomes simple in the limit of a very 
elongated condensate along z pl[ | . As shown in appendix 
0, the vortex 2D energy function £'^°(rj^;z) depends 
on r±_ only via the ratio r±/R±(z) where R±{z) is the 
local Thomas- Fermi radius. As a consequence El^'^°{Z) 
depends on Z through Z/Rz only, and its derivative is 
R±/Rz times smaller than the derivative of E^''S"^{Z). 
In the limit Rz/ R± tending to infinity, equation (|3q) re- 
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duces to 



dZ 



(39) 



Using the exphcit expression of E^^™ (Z) similar to ( |35| ) 
we obtain the condition 



(40) 



that is Z = Zc, the vortex starts bending at the elevation 
where the 2D critical rotation frequency equals the trap 
rotation frequency fJ. As a consequence the minimal ro- 
tation frequency to stabilize the bent vortex line is given 
in the lOz/ijJ± —^ limit by the 2D critical frequency in 
the central slice z = 0, see ( A68 ): 



fii = — i log 



,C+l/2 



rj 



(41) 



with 7] ~ 1.938 and C ~ 0.0884. This asymptotic predic- 
tion is plotted as a solid line in Figure || and is in good 
agreement both with the numerical minimization of the 
vortex energy functional and with the numerical mini- 
mization of the full Gross-Pitaevskii energy functional 

H- 

To conclude this subsection we point out two striking 
properties of the result (^l]). Firstly, it explains why in 
a cigar-shaped condensate, the bent vortex is first stabi- 
lized with an almost vanishing energy gap with respect 
to the 0-vortex configuration, at least much smaller than 
h'^Lu'^/fj,, see Figure U a,t fl = fli, the variational ansatz 
predicts an energy scaling as {h'^uj'j_/ ^)R±/Rz, see the 
expression for E^°''^^ with Zc/Rz — 0. This is very differ- 
ent from the 2D case, where the vortex, when first sta- 
bihzed, has a large and positive energy ~ h'^u;j_/ fi with 
respect to the vortex free configuration. Secondly, it is 
remarkable that the 2D stabilization frequency plays no 
role in the 3D case for the cigar-shaped condensate. This 
means that being a local minimum of energy in 2D slices 
does not imply that the vortex line is a local minimum 
of energy in 3D. As an example, it is possible in 3D to 
shorten the vertical segment of the vortex line, that is 
to introduce the bending at a slightly lower elevation, 
whereas this infinitesimal transformation has no equiva- 
lent in a purely 2D case. 



D. Saddle points of the simple energy functional 

Using a simple energy functional rather than the full 
Gross-Pitaevskii energy functional has allowed to under- 
stand the bending of the vortex line physically with the 
help of a simple piecewise variational ansatz. But much 
more can be done: as we now show, one can investigate 
not only the energy minimum but also possible saddle 
points of the energy functional. This will explain an in- 
triguing feature of the numerical results of Figure |l|: the 
angular momentum of the condensate has a discontinu- 
ous jump from the one-vortex branch to the no-vortex 
branch, which opens a gap in the allowed values of (Lz). 



The idea is to consider now a more general ansatz for 
the vortex line than the one of the previous subsection, in 
order to allow the vortex line to move towards the border 
of the condensate. We introduce the two parameter linear 
piecewise ansatz in the y — plane shown in Figure ^ 
a horizontal line linking {x = +oo,z ~ ^Z) to {x 
X, z = —Z), then a vertical segment of length 2Z parallel 
to z, and finally a horizontal line from (x — z = Z) 
to {x — +O0, z = Z). In other words, the vortex line has 
a C shape of adjustable distance X to the rotation axis 
and adjustable height 2Z. 

We consider the parameters of the ENS experiment, 
with A = uJz/u;± = 1/15. Eq.([2^) leads to a Thomas- 
Fermi chemical potential of /x = 7.621ioj±. We choose a 
rotation frequency = 0.448a;j^, which is a factor 1.2 
above the threshold value fli of (|4^). Figure |l^ shows a 
contour plot of the simple energy functional as function of 
the coordinates {X, Z) of the bending point of the vortex 
line. One recognizes first the global minimum of energy, 
marked with a M, with a position of the C shaped vortex 
line similar to the one of Figure |l^a: the vertical part 
of the vortex line is very close to the rotation axis, with 
X less than 0.05 times the transverse Thomas-Fermi ra- 
dius Second, one finds a saddle point in the energy, 
marked with a S, corresponding to a quite different posi- 
tion of the C shape: the vortex line is now far from the 
rotation axis, with X ^ 0.45i?x, and the half elevation 
of the C shape is of the order of R±^ . 

What is the physical meaning of this saddle point ? 
As it corresponds to an extremum of the energy func- 
tional (vanishing first order derivatives of the functional) , 
it represents a stationary shape for the vortex line in the 
rotating frame. In addition to this dynamic property, it 
has the following interesting energetic aspect: it corre- 
sponds to an energy minimum for a fixed value of the 
angular momentum per particle {Lz), as one can check 
explicitly for the two-parameter ansatz. 

These properties are not specific to the two-parameter 
ansatz. They apply for arbitrary shapes of the vortex 
line, but also for the exact Gross-Pitaevskii energy func- 
tional. Consider indeed the subspace of condensate wave- 
functions with a fixed angular momentum per particle lz ■ 
The Gross-Pitaevskii energy functional in the absence of 
rotational term —flLz is bounded from below (for g > 0) 
in this subspace so it has a minimum for some wavefunc- 
tion ipi^, with an energy E^^^^ilz). One can then show 
that the wavefunction ij^i^ is a stationary point for the full 
Gross-Pitaevskii energy functional, that is including the 
rotational energy term —VlLz, without imposing a fixed 
the angular momentum, provided that one takes the ro- 
tation frequency equal to = dE^^^ / dlz- One can show 
that it is however not necessarily a minimum of the full 
energy functional: it is a minimum if 



dVL 

dTz 



d^E 



dll 



> 



(42) 



otherwise it is a saddle point. 
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FIG. 11: For the two parameter linear piecewise ansatz for 
the vortex hne shape, contour plot of the simple energy func- 
tional as function of the coordinates {X, Z) of the bending 
point of the line. The trap aspect ratio is A = ujz/i^± ~ 1/15, 
leading to a Thomas- Fermi chemical potential /i — 7.62hu)±. 
The rotation frequency is f2 = 0.448ti;x ~ 1.2r2i. X and Z 
are expressed in units of the transverse Thomas-Fermi radius 
R± of the condensate. The contour plot is restricted to a 
bending point inside the Thomas-Fermi profile of the conden- 
sate, -I- A^Z^ < The global minimum is marked with 
the letter M, and the saddle point is marked with the letter 
S. The elliptic isocontours close to the Thomas-Fermi border 
correspond to a maximum of energy. 



We can exemplify these properties with our simple two- 
parameter ansatz. We have plotted in Figure 12 the mean 
angular momentum per particle of the energy minimum 
and of the saddle point as function of the rotation fre- 
quency f2. One finds that the branch of minimum has an 
angular momentum ranging from 0.42/i to ?i, whereas the 
branch of saddle point fills the gap of angular momentum 
from to 0A2h. One also finds that the angular momen- 
tum is an increasing function of the rotation frequency 
on the branch of minimum, whereas it is a decreasing 
function of Q on the saddle point branch, in accordance 
with the general condition (p^. 

In an experiment where the angular momentum of the 
condensate is gradually decreased, starting from a value 
close to h, the vortex first follows the minimum energy 
branch: the distance of the vortex line to the rotation 
axis X remains much smaller than Rj_, whereas the vor- 
tex half-height Z and the rotation frequency 57 of the line 
in the lab frame decrease with time. Then the angular 
momentum reaches a critical value: the minimum in the 
{X, Z) parameter space merges with the saddle point. 
When the angular momentum is further decreased, the 
vortex line follows the saddle point branch: the distance 
X increases and becomes of the order of R± , as the vor- 
tex moves towards the boundary of the condensate; Z 
takes values as small as R±, in which case the validity 
of our simple analytical approach becomes marginal, and 
the rotation frequency of the vortex line increases. This 
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FIG. 12: For the two-parameter ansatz and the simple en- 
ergy functional for the bent vortex line, mean angular mo- 
mentum per particle for the energy minimum (solid line) and 
the energy saddle point (dashed line) as function of the ro- 
tation frequency fl. The trap aspect ratio is A = 1/15 and 
the chemical potential is jj. — 7.62hiL}±. The points marked 
with M and S correspond to the rotation frequency of Figure 
|ll| . The vertical line corresponds to the analytical prediction 
for Sli. The rotation frequencies are in units of uj± and the 
angular momentum is in units of h. 



scenario was very recently observed at ENS ||l8| . 

To conclude this subsection, we can emphasize again 
the analogy of the 3D bent vortex problem and the 2D 
single vortex problem. In 2D, the minimum of energy 
in the rotating frame with no constrains! on the angular 
momentum corresponds to the vortex at the trap center, 
with an angular momentum h, or to the vortex at infinity 
with a vanishing angular momentum. If one minimizes 
the energy for a fixed angular momentum Iz, one finds 
an off-center vortex core, shifting from the trap center 
to infinity as Iz is varied from h downto zero: these off- 
center solutions do not satisfy the criterion (^) and 
are very analogous to our saddle point solutions. 



V. FINITE TEMPERATURE FLUCTUATIONS 
OF AN OTHERWISE STRAIGHT VORTEX LINE 

A. Method: Bogoliubov approach around the 
straight vortex steady state 

In order to evaluate the contribution of thermal fluc- 
tuations on the contrast, we study the limiting case of a 
straight vortex. For the ENS trap parameters, this means 
that 1/A < 7 and in all the following, we illustrate our 
discussion using the case (A = 1/5, = 0.5oj±). In this 
regime, the system has a rotational symmetry around the 
z axis so that the numerical problem is effectively 2D in 
cylindrical coordinates. We use the ?7(l)-symmetry pre- 
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serving Bogoliubov approach described in [g3|. In this 
way, the problem of the spurious mode of the condensate 
is avoided. The Bose field is expanded as 



0.04 



7^(f ) = (t){r)a^ 



A^Y.^bkUk{f)+blvl{r). (43) 



(j) is the condensate wave function normalized to unity. 
The operator is the annihilation operator in the con- 
densate mode and the almost unitary operator — 
" l)~^^'^a^ gives the phase factor of the field in the 
iq]. The modal functions Uk,Vf^ 



condensate mode 
are normalized like 



r \Uk\ 



\vk\ = 1. They are ob- 



tained from the usual modal functions of the Bogoliubov- 
de Gennes equations after orthogonalization of u^. and of 
Vk with respect to (p. The Bose operators bk annihilate 
one quasi-particle in the mode k but conserve the total 
number of particles. The index k = {n,l,s} denotes the 
quantum numbers of the mode linked to the symmetry 
of the system: s = 1 (or —1) for symmetric (or antisym- 
metric) modes with respect to the plane z — 0, Ih is the 
angular momentum with respect to the condensate and 
the integer n is the radial quantum number. 

Concerning the transverse direction x — y, we expand 

(f>, Uj^, Vk on the harmonic oscillator basis |3>[|^ | of pul- 
sation uj±_ {m'h is the angular momentum and n', the 
radial quantum number). For example, in this basis we 
have 



(f) 



ri'=0 



(44) 



Numerically, the basis is truncated: the spatial grid 
along z is surrounded by infinite walls, and also the num- 
ber of wave functions is limited (for the value 
A = 1/5, the harmonic oscillator basis in our compu- 
tation contains all the wave functions of energy less than 
42 fiLo±). The choice of the grid and the discretization of 
the Laplacian along z has been made so that the first 160 
energy levels e„ of the pure ID harmonic oscillator are 
recovered with an error |e„ — e°'^'''^*| < 10^^. In this part 
we have computed the condensate wave function using 
an imaginary time method and the convergence criterion 
® (typical values of 7 are of the order of 10^^'^). 



B. Expectation values of some observables 

The mean density of atoms out of the condensate is 
obtained straightforwardly with the usual expression 

Poxc(r ) = + ) + \vk\^{r )) , 

k k 

(45) 

where 
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FIG. 13: Mean column density (n) associated to the param- 
eters (A = 1/5, f2 = 0.5ti;j_), as function of the distance r 
to the rotation axis. The temperature T ranges from to 
0.4rc in steps of 0.05 Tc, from bottom to top. r is in units 
of the harmonic oscillator length ax ~ {h/mui^^^^^ and (n) 
is in units of Na^ — NmLUj_/h. The number of particles is 
= 1.4 X IC and the chemical potential is ^ ~ 9.87?iaJx. 
The inset is a magnification. 



is the Bose occupation factor and efc is the energy of 
the Bogoliubov mode k in the rotating frame. Because 
of numerical constraints we have limited the sum over 
states of energy less than 20 1ilj±^ in the rotating frame, 
and we have studied configurations with temperatures 
less than 11.4 hu)±/kB that is less than 0.4Tc, where 

1/3 

/cbTc = h (a;2W^7V/C(3)) is the ideal Bose gas critical 
temperature. For comparison with experimental results, 
we extract the mean column density (n) in the x — y 
plane: 

{n){x,y) = Jdz [(TV - (7Vexc))|0P(r) -f Poxc(r)] , 

(47) 

where (A^cxc) = j d'^rpc^d^^) is the mean number of 
atoms out of the condensate. We plot this column den- 
sity as function of the distance to the z axis for different 
values of the temperature in Figure |l^. At zero temper- 
ature the effect of quantum depletion at the center of the 
vortex line is not observable on the scale of the figure. 
Hence quantum fluctuations of the vortex line are clearly 
not responsible for the low contrast measured in ENS 
experiments. 

To be more quantitative, we denote (n)™'" the value 
of the mean column density at the center of the trap and 
^^(jmax .j-j^g maximal value of the mean column density. 
Then, the mean contrast is defined as 



C = 1 



(48) 



The modal functions are of the form 



Uk = l/(exp(efc/fc_BT) - 1) 



(46) 



Uk{r) = ?7fc(r, z) exp + 1) 



(49) 
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FIG. 14: Dependence of the mean contrast C — 1 — 
(n)""°/{n)™'"' with the temperature for the same parameters 
as in Figure |l^. 

Vk{r) = Vk{r,z)ex_v[i{l-l)e] , (50) 

If / + 1 7^ the function Uk vanishes in r = because of 
the centrifugal barrier. For the same reason Vk{r = 0, z) 
vanishes for I — 1 0. As a consequence only modes 
with / — ±1 contribute to (n)""". Furthermore, in the 
rotating frame the energy of the quasi-particles is given 
by Eq.(^ and the most populated modes are the lowest 
energy kelvons, characterized by / — —1, which have a 
negative energy in the absence of rotation 

111- In 

the example presented in this section, the lowest kelvon 
mode localized at the extremities of the condensate has 
an energy in the rotating frame given by i? ~ 0.018 hLu±. 
Figure |l^ shows that the dependence of C as a function 
of temperature is almost linear. This behavior can be 
understood from the fact that the lowest energy modes 
have a semi-classical character {ek ^ ksT)^ with an 
occupation number depending linearly on temperature 
rife ~ kBTlek. 

Figure |lj shows that temperature contributes signif- 
icantly to the observed contrasts. In particular, for a 
temperature of the order of the chemical potential, which 
is a typical situation in experiments, we find that the 
contrast is 65%, not far from the contrast value of ^ 
50% observed at ENS. Note that some zero temperature 
bending of the vortex line is expected to occur for the 
parameters of ENS, see Figure which adds up to the 
thermal fluctuations and further reduce the contrast in 
the experiment. 



C. Mimicking a single experimental run: a 
Glauber P method 

Let us insist now on the fact that expectations values 
calculated in the previous section are not sufficient for 



a quantitative understanding of experiments. They just 
give an order of magnitude of what is observed. Indeed, 
for example, for the particular value of A that we have 
chosen, the straight vortex is at the edge of the existence 
domain and one expects large fluctuations of the phase 
and of the density due to the soft core mode and also to 
the emergence of low energy surface modes. As a conse- 
quence, observables obtained from a single measurement 
of the system may differ notably from their mean value. 
Hence, the aim of this section is to evaluate the fluctua- 
tions induced by the soft modes. Calculations of the pre- 
vious subsection have shown that quantum fluctuations 
do not contribute to the contrast at finite temperature, so 
that semi-classical field approaches are good tools to an- 
swer this problem. We use for that purpose the Glauber 
P method of quantum optics . 

This approach is made simple by the fact that 
the A^-body density operator D is the exponential of 
^H^ogl {ksT), where the Bogoliubov Hamiltonian is a 
sum of decoupled harmonic oscillator terms, iJeog — 
Eq + J2k^kbl.bk- We introduce the coherent state of 
quasi-particles \{Pk}) = \(3ki, fik2, Pks ■ ■ ■) , where \f3k) is 
the normalized eigenstate of the annihilation operators 
fefe with the eigenvalue Pk ■ Then one defines the Glauber 
P representation of D: 

i>= /n ''''^'''"^' ^({/^^'})i{/^w)({/3.4i. (51) 

k 

The explicit calculation of P is possible, as D is Gaus- 
sian in the 6fe's this gives the following Gaussian 
distribution: 

pm'})^ll-^-A-^) ' (52) 

rik \ Uk J 

where nk is the mean number of quasi-particles in the 
mode k given in Eq. (^). This Gaussian distribution is 
easily sampled. As P is positive one can imagine that a 
given experimental realization of the gas is in the state 
|{/3fe}) where the complex numbers /3fe vary randomly 
from one realization to the other. In the following, we 
want to determine the density and the velocity field of 
the atomic gas for a given Monte Carlo realization of the 
{Pk}. 

The A^-body distribution function corresponding to a 
single term of the statistical mixture ( pT| ) 

p(ri,...,rW) = m'}\i^\ri)..4\rN) 

x^{TN)...i^{ri)\{(3k'}) (53) 

is not easy to calculate for the interacting Bose gas as 
ip is a. superposition of bk and b\ so that the product of 
field operators in ( p3| ) is not normally ordered in terms 
of the bk- We perform the following approximation 

bim'})^mPk'}). (54) 

The error has a root mean square norm equal to unity, 
hence the approximation is good for modes with a large 
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z(m) 



FIG. 15: Representation of the x coordinate (solid line) and 
the y coordinate (dashed line) of the vortex line as function 
of the position z on the rotation axis, for a siiigle stochastic 
realization. The parameters are as in Figure llSl with T = 
0.2Tc. Note that the transverse Thomas-Fermi radius is of 
the order of 6/im. 



occupation number, bad for empty modes. In this way 
we describe correctly the fluctuations due to finite tem- 
perature, but not the quantum fluctuations existing even 
at zero temperature. As we have shown previously, this 
is fine in the dilute limit {pa^Y^"^ <C 1, where quantum 
depletion is small. The approximation amounts to taking 



i,{f)\{(3k'})^^{r)\{Pk'}) 



(55) 



with a classical field 



>^ 




FIG. 16: For an individual Monte Carlo realization of the 
Bose field tp, position of the vortices in the velocity field of 
the gas in the plane z — 0, for a temperature T — 0.2Tc. The 
vortices with a positive charge are indicated with a cross, the 
vortices with a negative charge are indicated with a triangle. 
The circle indicates the Thomas-Fermi border of the conden- 
sate. The X and y coordinates are in units of {Ti/muj_i_Y^^ . As 
a consequence of thermal excitation of the low energy surface 
modes, satellite vortices appear at the border of the conden- 
sate, around the central vortex. The physical parameters are 
as in Figure 



where iVo is such that the norm squared of is equal 
to the total number of particles N In this approxi- 
mation, a single realization is now in the coherent state 
|{'0}) for the field ■0. In this case the TV-body distri- 
bution function of the atoms for a single experimental 
realization is factorized: 



(57) 



In this paragraph, we describe temperature effects in 
a single stochastic realization of "0 at fcsT — 0.2kBTc ^ 
^fj,. First, we have extracted the shape of the vortex line: 
Figure 15 shows the x and y coordinates of the vortex line 
as function of z. The soft core modes localized at the two 
extremities |l^ of the condensate are thermally ex- 
cited and as a result push the core away from the rotation 
axis. This is a temperature induced bending. Second, for 
the same stochastic realization. Figure |l^ represents the 
location of the vortices in the velocity field of the gas in 
the plane z = 0. It shows that satellite vortices, which are 
not easily observable in the density profile, appear at the 
border of the condensate. The presence of these vortices 
is due to the excitation of surface modes of low energy 



and relatively high angular momentum (recall that for 
A — 1/5, il — 0.5llj±, the system is at the border of the 
thermodynamic stability domain). 

We address now the problem of the integrated den- 
sity contrast for a single stochastic realization. In the 
previous single stochastic realization at a temperature 
T = 0.2Tc, the value of the contrast is around 90% which 
does not coincide with the mean contrast C ~ 85% de- 
rived in the previous subsection. This disagreement sug- 
gests strong fluctuations of density along the z— axis. To 
confirm this idea we have computed 10000 stochastic re- 
alizations of the column density at the center of the trap. 
We have reported in Figure |^ the corresponding his- 
togram. This figure shows effectively that the probabil- 
ity law is far from being Gaussian: the long tail on the 
right side is due to ID character of the excitation modes 
of the vortex line, the so-called kelvons, similar indeed 
to the fluctuations of the number of condensate parti- 
cles in a ID Bose gas This is another indication of 
the importance of fluctuations in the properties of this 
system. 

To be complete we now consider the effect on the con- 
trast of the ballistic expansion of the cloud performed 
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FIG. 17: Histogram associated to the single realization 
column density n at the center x — 0,y — 0, for 10^ re- 
alizations. This figure shows that the probability law of 
this observable is strongly non Gaussian, n is in units of 
N/a]_ = Nmuj±/h. The same parameters as in Figure ^ 
were used, with T = 0.2Tc. We have observed that this non- 
Gaussian character subsists even at a temperature as low as 




-5 5 

y (units of a ^, rescaled) 

FIG. 18: Comparison of the column density for a straight 
vortex (A — 1/5) at finite temperature (T = 5.71 Tiuj±/kB) 
before and after the release from the trap (time of flight ~ 
30ms), as a function of the rescaled y coordinate. The same 
parameters as in Figure ^ were used. 



in the experiment. As in subsection HID we integrate 



numerically the rescaled Gross-Pitaevskii equation, with 
a classical field ijj including thermal fluctuations. As we 
see in Figure ^ the rescaling now absorbs to a lower ex- 
tent the effect of the ballistic expansion, but the density 
contrast is only weakly changed. 



VI. CONCLUSION 

We have understood why the vortex line can bend in 
a steady state cigar-shaped condensate rotated at fre- 
quency f2 around its long axis z. In the Thomas-Fermi 
regime the cigar-shaped condensate can be viewed as a 
collection of slices parallel to the x — y plane, each slice 
corresponding formally to a 2D rotating condensate. For 
each 2D condensate one defines as usual the critical rota- 
tion frequency Vl^ above which it is energetically more 
favorable to have the vortex core at the trap center rather 
than at infinity. As the 3D condensate density is inho- 
mogeneous along z axis due to the harmonic confinement 
the local f2^^ is minimal in the plane z = and maximal 
at the end points of the cigar. 

The vortex line then uses the following strategy to min- 
imize its energy: it follows the rotation axis z in the ele- 
vation interval where the trap rotation frequency is larger 
than the local fi^^'s, and moves away from the rotation 
axis to infinity where O becomes smaller than the local 

This leads to the analytical prediction that the minimal 
rotation frequency f2i required to stabilize the bent vor- 
tex in a cigar-shaped condensate is equal to ^^{z — 0), 
that is the 2D critical rotation frequency corresponding 
to the central slice z = 0. 

Another analytical prediction is that, for > Jli, the 
energy presents a saddle point, in addition to the previ- 
ously discussed minimum. This saddle point is a station- 
ary configuration of the vortex line in the rotating frame, 
with an angular momentum smaller than the one of the 
energy minimum. The vortex line in the saddle point has 
also a different position with respect to z axis than in the 
minimum energy configuration: the vortex line, rather 
than being on axis, is shifted away from the rotation axis 
by an amount of the order of the transverse Thomas- 
Fermi radius R±^ , and it has a much smaller length along 
z, also on the order of This description matches very 
recent experimental observations of bent vortices at ENS 

To test the analytical predictions for the energy min- 
imum we have performed a full numerical minimization 
of the Gross-Pitaevskii energy functional with the effi- 
cient conjugate gradient method for typical parameters 
of the ENS experiment. We have constructed in this way 
a phase diagram giving the existence domain of a single 
vortex configuration as a local minimum of energy, as 
function of the trapping potential aspect ratio and of the 
rotation frequency. The numerics confirm the analyti- 
cal prediction for fJi, and give the rotation frequency ^2 
above which surface modes of the condensate are desta- 
bilized and several vortices enter the condensate. 

We have also studied the effect of thermal fluctuations 
of an otherwise straight vortex line in a cigar-shaped con- 
densate. We find values of the mean density contrast of 
the vortex hole close to the experimental results of [Q. 
Using the Glauber P representation of the density oper- 
ator of the gas for the field of Bogoliubov quasi-particles. 
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we can predict what happens in a single reaUzation of 
the experiment: the vortex hne experiences some bend- 
ing, due to the thermal population of low energy kelvon 
modes localized close to the end points of the condensate. 
The ID nature of the kelvon modes leads to a remarkable 
non Gaussian character of the vortex line fluctuations. 



then finds an energy difference between the vortex free 
configuration (of energy Eq) and a configuration with a 
vortex of the order of 



E ~ Eq 



2, ,2 



(A2) 
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APPENDIX A: SIMPLE DERIVATION OF A 
VORTEX LINE ENERGY FUNCTIONAL 

We present here a derivation of an approximate energy 
functional for a single vortex line in a Thomas-Fermi con- 
densate subject to a harmonic potential rotating around 
one of its eigenaxes. Such a derivation is available in 
the literature |l^ but it is rather complex as it directly 
includes the effect of curvature of the vortex line in the 
energy functional, and it applies for an arbitrary aspect 
ratio of the condensate. The derivation that we propose 
has the advantage of simplicity as (i) it is specialized to 
the case of a very elongated condensate along the rota- 
tion axis z in a cylindrically symmetric trap, that is the 
atomic oscillation frequency along z is much smaller 
than the oscillation frequency uj± in the transverse xy 
plane, and (ii) it immediately supposes that the explicit 
dependence of the energy functional on the curvature 
of the steady state bent vortex line can be neglected. 
Point (ii) allows to assume locally that the vortex line 
is straight, which greatly simplifies the derivation of the 
energy functional. By comparison with the more general 
derivation of [0 we have checked the validity of point 
(ii) (2|. 

Furthermore we assume in our derivation that the ro- 
tation frequency Q is of the order of 



A* 



(Al) 



where lo^_ is the oscillation frequency of the atoms in the 
xy plane and /i is the chemical potential. It is indeed 
in this range of rotation frequencies that the single vor- 
tex configuration is first stabilized, and this also greatly 
simplifies the derivation of the energy functional. One 



The first step in the derivation of the approximate en- 
ergy functional is to rewrite the full Gross-Pitaevskii en- 
ergy functional (|l|) in terms of the modulus and the phase 
of the condensate wavefunction 0(r ): 



(f) 



,1/2 



(f)e' 



iS{r) 



(A3) 



where n is the probability density normalized to unity. 
At this stage it is convenient to introduce the so-called 
velocity field of the condensate: 



vir) = — gradS'(r). 
m 



(A4) 



Note that v corresponds to the velocity field in the lab 
frame, the velocity field in the rotating frame being 
V ~ Q. f\r. Inserting (A3) into the Gross-Pitaevskii en- 



ergy functional (|l|) leads to an expression that we split 
for convenience in two contributions, one linked to the 
modulus and the other one linked to the phase: 



E — E-^nod ^ -^phaso 



(A5) 



= d 



ig-axd^/^f + Un+ -NgnHPS) 

2m 2 ■ ^ 



E, 



phase 



—mv"^ — ^l ■ (r A mv] 
2 ^ ' 



(A7) 



As already mentioned, the trapping potential U is axi- 
symmetric with respect to z: 

U=^m[Ljl{x^+y^)+u;lz^]. (A8) 

We also give the conditions on n and v ensuring that is a 
local extremum of the Gross-Pitaevskii energy functional: 



= div 
1 

fi — —m 



n{v — r2 A r ) 
+ f/ -I- Ngn 



2m Jn 



(A9) 

mfi f\r ■ V 

(AlO) 



where [i is the chemical potential. These conditions are 
simply the time independent Gross-Pitaevskii equation 
written in the modulus-phase representation. 



1. Energy from the modulus 

In the present Thomas-Fermi regime the contribution 
depending only on the modulus can be evaluated along 
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the lines of our previous work ||]. One splits the den- 
sity in a slowly varying envelope rigiow and a function 
representing the density hole due to the vortex core: 



n{f) = nsiowir)f^{r)- 



(All) 



The envelope ngiow varies at the scale of the trans- 
verse Thomas-Fermi radius of the condensate B.± = 
whereas / varies at the scale of the di- 
ameter of a vortex core, that is the healing length ^ = 

We now determine the function / from the requirement 
that it deviates significantly from unity only at a distance 
at most a few ^'s from the vortex core. As such a length 
scale we neglect the spatial variation of U and of rigiow; 
we further check that in the range (Al) of rotation fre- 
quencies, the rotational velocity term 17 A r is negligible 
as compared to the vortex velocity field v, which allows 
to negl ect th e A r terms close to the vortex core both 
in (|M and (|aTo|): 



\nAr\ flRj 



n 



< 1. 



(A12) 



Furthermore the minimal radius of curvature of the vor- 
tex line is found in the subsequent calculations to be of 
the order of R± much larger than ^. After all these sim- 
plifications the function / is found locally to solve the 
well-known Gross-Pitaevskii equation for a straight vor- 
tex line in an infinite spatially homogeneous condensate 
1^ , this fictitious homogeneous condensate having a par- 
ticle density given by Nusiow evaluated on the vortex line. 
We therefore take for /: 



/(f) = F(d/6oc) 



(A13) 



where d is the distance of f to the vortex line and 



s,loc 



y/mgNusiowiro) 



(A14) 



is the local healing length at the position of the vortex 
line po| . The function F{u) does not depend on any 
physical parameter. It can be obtained from a numerical 
solution of the reduced Gross-Pitaevskii equation for a 
vortex in a homogeneous condensate (see equation (2.84) 
of j|]). In the large u limit, its deviation from unity tends 
to zero as 0(1/m^). For moderate values of u it is well 
approximated by 



F{u) ~ tanh(0.7687M). 



(A15) 



The slowly varying envelope Ugiow is obtained from 
by removing the short range mf^/2 term already 



included in / and by neglecting the quantum pressure 
term cx Ay^/y^ in the spirit of the Thomas- Fermi ap- 
proximation. We write Ugiow as utf + Sn where 



JiTFir) 



/^o - U{r) 
Ng 



(A16) 



is the usual Thomas-Fermi approximation for the proba- 
bility density in the absence of vortex and 



NgSn = Sfi + Q A f ■ mv. 



(A17) 



The density correction 5n includes the rotational term 
and the deviation J/i of /i from the vortex free Thomas- 
Fermi chemical potential /ip. Both contributions are of 
the same small order. Using the estimate v h/mR±^ 
we find that the rotational term leads to 



nn 



Mo 



(A18) 



We can also estimate by multiplying (AlO) by n and 
integrating over the whole space: we find a contribution 
to 5^ involving the rotational term. 



(A19) 



of the order of Ml. A second contribution comes from the 
fact that the vortex line digs an empty tube of volume 
^ Rz^^ in the condensate of volume ~ R\Rz, where Rz 
is the Thomas-Fermi radius along z, so that fi has to 
differ from /iq by a relative amount 



Mo 



R\Rz 



Mo 



(A20) 



to ensure that n is normalized to unity. 

We now proce ed with the calculation of i?mod , inserting 
the ansatz (All) in ( A_6). We calculate first the harmonic 
plus interaction potential energy part of i^modj then the 
kinetic energy part of i^mod- 

In the harmonic plus interaction potential energy 
terms, we use the identity — [p — 1) -I- 1 and we 
collect the terms in powers of — 1: 



pot 

(0) ^ 

pot 

p(i) _ 



pot 



(1) 



E, 



(2) 



E, 



pot ~ pot 



low 



(A21) 
(A22) 



E, 



pot 

(2) 
pot 



dV(A^gnsiow + U)n,xoAf - 1) (A23) 



(A24) 



The zeroth degree term is of the order of magnitude of 
/i so that one has to include the deviation of nsiow_from 



nxF to first order in order to get the leading term ( A2 ) in 
the vortex energy. Using the fact that NgriTF + = Mo 
we obtain 



Epot ~ E' 



TF 



Mo 



£fSn + o{n^ujl/^i), (A25) 



where we have introduced the Thomas-Fermi approxima- 
tion to the energy Eq of the vortex free configuration: 



En 



Etf = / d^r [ ^nlp 



U riTF 



(A26) 
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The first degree term E^^^ and the second degree term 
E^J^ are of the order of Ti^w^/^o as |/^ — 1| is close 
to unity in a cyHnder of volume Rz^^ and is negligible 
outside. We can therefore approximate risiow by nxp in 



E^ll and eI'I 



We then get for E' 



(1). 

pot ■ 



E, 



(1) 

pot 



Jd^rfiTFif-l) 



= -Ho J d^rSnf ~ -^0 j d?r5n (A27) 

where we have used the normalization of rixF and of 
^siow/^ ~ (?^TF + 5n)p to unity. We thus see that at 
the present order of the calculation E^^^ compensates the 
term linear in 5n in -Epoj. We finally get: 



E, 



pot 



(A28) 



In this integral we introduce a local system of cylindrical 
coordinates (p, 9, Z) with a vertical axis Z tangent to the 
vortex line, see Figure [l9|, the first coordinate p being the 
distance to the vortex line. We then approximate nxp by 
its local value on the vortex line. The angular integral 
over 9 gives a factor 2tt. The function — 1 depends 



on the distance p to the vortex line, see equation ( Al3| ), 
and the corresponding integral over p can be extended 
to infinity as (/^ — 1)^ tends rapidly to zero far from 
the vortex line. We parametrize the vortex line by its 
curvilinear abscissa s. We then realize that in the integral 
over Z we can make the reinterpretation dZ = ds so that 



Epot — Etf 



(A29) 



J ds ^n2p(fo(s)) ^ 27Tpdp [f2(p/6oc(s)) - 1]' 

where ro(s) is the position of the vortex line at abscissa 
s. Finally rescaling p by 6ocin the integral over p and 
replacing risiow by titf in ( A14) leads to 



-Epot - Etf ~ Aq ds nTF(fo(s)) 

J m 



where the constant factor Aq is 



Aa= udu [F{uf ~ l\ 



(A30) 



(A31) 



In the kinetic energy term of -Emodj we neglect the 
spatial derivative of the slowly varying envelope Ugiow 
and of the local healing length ^, as they both vary on 
a length scale R± . We calculate the gradient of / in the 
local system of cylindrical coordinates (p, 9, Z) of Figure 

m 



(grad/)^ 



1 

Slot 



F'ip/^loc). 



(A32) 



This function vanishes at p ^ ^joc from the vortex line so 
that we can replace Ugiow by its value on the vortex line 



z /N 
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FIG. 19; Local frame with Cartesian coordinates X, Y, Z 
around a point fo{s) of the locally straight vortex line. Local 
axis Z is tangent to the vortex line and has the orientation of 
the local vorticity, so that the unit vector along Z coincides 
with the vector t{s) tangent to the vortex line defined in the 
text. X, Y are arbitrary Cartesian coordinates in the plane 
orthogonal to Z. One then defines the corresponding cylin- 
drical coordinates p,9,Z of an arbitrary point M, e.g. p is 
the distance of AI to the vortex line. 



and extend the integral over p to infinity. The integral 
over the angle 9 gives a factor 27r. As in the previous 
paragraph we use dZ = ds where s is the curvilinear 
abscissa on the vortex line, we rescale p by ^loc in the 
integral over p and we replace rigiow by the vortex free 
Thomas- Fermi expression utf- We then obtain for the 
kinetic energy part of Emoci: 



Tjikin 
-'^mod 



Ai 



ds^^riTFirois)) 
m 



where the constant factor Ai is given by 



Ai - 



udu F ■^(u). 



(A33) 



(A34) 



To summarize the vortex energy stored in the modulus 
of the condensate wavefunction is approximated by 



^mod 



E' 



ttTl 



TF - {Aq + Ai) / ds nTF(ro(s)). (A35) 

m 



We estimate E^mod — Etf to be indeed of the order of 
magnitude of h^LU'^/po, as expected from (A2), by tak- 
ing a length for the part of the vortex line inside the 
condensate and a typical value for the Thomas- Fermi en- 
velope Utf ~ ^/RzR\- 
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2. Energy involving the velocity field 

The energy term Sphaso is quite difficult to evaluate in 
its present form (A7): since the velocity field tends to 
zero slowly away from the vortex line (as H/md where 
d is the distance to the vortex line) , all the parts of the 
Thomas-Fermi volume of the condensate give approxi- 
mately the same contribution. A direct 3D integration 
cannot be performed as the velocity field far from the 
vortex line is not known explicitly ISll , with the notable 
exception of a straight vortex line on axis z. 

The trick that we borrow from px| ] is to transform the 
volume integral ( A7) in a line integral, by 'integrating by 
part' and using the fact that curli; is known exactly: 

curlw(f) = — / dst{s)S{f- fo{s)) (A36) 
m J 

where s is the curvilinear abscissa along the vortex line 
and t is the unit vector tangent to the vortex line and 
oriented in the direction of the local vorticity {t has a 
positive component along z for a positive rotation fre- 
quency f2). Physically one can use an electromagnetic 
analogy: curl?7 corresponds to a linear current, that is 
a perfectly filiform charge current, of intensity 27r?i/rn 
circulating in the vortex line. 

The 'integration by part' to be performed relies on the 
following vectorial identity 



div(A AB) = B ■ curl A - A ■ curlB. 



(A37) 



Integrating this identity over whole space and using Os- 
trogradski's formula we arrive at the desired identity 



d-^fA ■ curli? 



d-^fB ■ curl A 



(A38) 



where we have assumed that the total flux of the vector 
A A B vanishe s thr ough a surface at infinity. 

We apply ( A38 ) first to transform the contribution 



E^tsc of«' to (^: 



kin 

phase 



d f—mnv ' 
2 



We have to choose A = v so that curlw appears in the 
left hand side of (|A38|). We then set flsl: 



curls 



(A40) 



This is an approximation as nv does not have strictly 
speaking a vanishing divergence. The vector which has 
an exactly vanishing divergence is the probability current 
in the rotating frame, n{v — ^Af), see (|A9|). Fort unately 
we can repeat the reasoning of subsection (Al). At a 



distance of up to a few healing lengths ^ from the vortex 
core, the rotation term J7 A f is negligible as compared 
to the velocity field v, see (A12) so that (A£) reduces 
to div(n'y ) ~ 0. At a distance much larger than ^ from 



the vortex core the density n can be approximate d by 
the vortex free Thomas- Fermi density tt-xf so that (A9) 
reduces to 



nTF{v — n Ar) 



0. 



(A41) 



As the Thomas-Fermi density nxF is rotationally sym- 
metric with respect to the rotation axis z, we have p3] 



div 



riTF^ A f 







(A42) 



so that (A41) reduces to div{nv) ~ 0. This finally leads 
to the desired line integral reformulation for iS'^J" ,„: 

c piiase 



TT^kin 
phase 



wh / dsB{ro{s)) ■ t{s). 



(A43) 



The calculation of B remains a challenge. Formally B 
can be considered as the static magnetic field created by 
a current proportional to nv so that we have the Biot 
and Savart formula [tM|: 



B{f) 



1 

47r 



,3 ^, nv (r ') A (r — r ') 



d^f 



?/|3 



(A44) 



but this requires in principle the knowledge of the veloc- 
ity field V everywhere ||3^. In practice the problem is 
simplified by the fact that one needs to know B on the 
vortex line only, f = f'o (s) fo r any fixed s, and by the 
fact that the integrand (A44) tends rapidly to zero for 
increasing |ro(s) — so that it is sufficient to know the 
velocity field v close to the vortex line. As mentioned in 
the beginning of this appendix we assume that the vortex 
line is locally straight around the point ro(s). We then 
approximate v by the velocity field of a straight vortex 
in a homogeneous medium. In the local system of cylin- 
drical coordinates (p, 9, Z) defined in Figure |l^ we thus 
write 



v{r') 



m p 



(A45) 



(A39) In this local frame the vector r' — ro(s) is equal to pep - 



Zez and the unit vector ez tangent to the vortex line 
actually coincides with t (s) so that one has 



t».[t/(f')A(fo(s)-f')]=^-. 

TO 

This leads to the rather explicit expression 

n(r ') 



kin 
phase 



4to 



ds / 



\ro{s) 



?/|3 ■ 



(A46) 



(A47) 



To calculate (A47) we write n = risiow/^ — «tf(/^ — 
1) -I- nxF as in subsection Al. This leads to a splitting 
of -EpJ'asc ill two pieces. 

The piece involving — 1 is re-expressed in terms 
of the local cylindrical coordinates (p, 9, Z) of the local 
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XYZ frame of Figure n9^: 



^phasG(I) 



4m 



4m 



^3^, nTF{r')[f{r')-l] 



ds / d^R 



,p"TF(r')[^'(p/eioc(s))-l] 



(p2 + Z2 + £2)3/2 



+ 00 



dZ 



1 



(Z2 + p2 + £2)3/2 p2 + , 



(A48) 



^phasG(I) 



dsnTF(n)(s)) 



(A49) 

A more explicit form will be given in the next subsection. 

The remaining piece of -Eph^gp involves the Thomas- 
Fermi envelope only: 



-^phasc(II) 



4m 



ds 



TiTFir') 



(|ro(s)-r'|2 + e2) 



3/2 ■ 



dz' 



1 
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the integration over z' to infinity and we use ( A48| ) to get 



where we introduced an arbitrarily small e in the de- 
nominator to prevent a divergence of the integral. The 
integrand as a function of Z is tending to zero as soon as 
\Z\ exceeds a few times p. As p is limited to a few times 
Cioc(s) by the function _F2 — 1, the integrand becomes neg- 
ligible as soon as \Z\ exceeds a few times £,\oc{s), which 
allows to approximate the slowly varying Thomas-Fermi 
envelope titf by its value on the vortex line and to extend 
the integration over Z to infinity. Using 



and extending in the resulting integral the integration 
over p to cx) we obtain the result 



• 



(A50) 

We simplify this expression by taking advantage of the 
cigar shaped nature of the condensate. We first integrate 
(A47) over z' . As in the previous paragraph we use the 
fact that the integral 



converges to its C — +oo value as soon as C exceeds a 
few times |fb_L(s) — r'^\, where rj_ is the projection of 
the vector r in the xy plane. The range of |r_L| is limited 
to the transverse Thomas-Fermi radius of the condensate 
by the presence of nTF(?^') in the integrand. The range 
of |ro^(s)| exceeds R±_ for a bent vortex line as the vortex 
line gets out of the Thomas-Fermi profile of the conden- 
sate; however the contribution to the energy of the vortex 
line segments at a di stan ce exceeding a few i?x 's becomes 
much smaller than (Al) and is hence negligible, see the 
appendix We can therefore assume that |ro-L(s) — f'^\ 
is at most a few times in (A50) so that the inte- 



gral over z' converges over a distance of R±. At such 
a length scale along the rotation axis z, titfI^;', y', z') is 
almost constant for a cigar shaped condensate and can 
be approximated by riT-pix' ..y' zq{s)). We then extend 



kin 

phase 



(11) 



2m 



ds I dx' dy' 



nTF{x',y',zo{s)) 



\ra±{s)-f'\^ + e^' 

(A52) 

The resulting integral over x' , y' can be calculated exactly 
and was already encountered in the 2D calculation of [|j . 
We will give the result in the next subsection. 

Fin ally we apply the 'integration by part' technique 
(|A3^ ) to the last term of -Bphaso, the rotational energy 
term: 



phase 



^ d 



— • r A mv 



d^fv-{nrif\f). 



(A53) 

Note that this term is simply —il{Lz), where (Lz) is the 
angular momentum per particle along z. Calculation of 
this term is considerably simpler than E^^^^^. First one 
can approximate n by the Thomas-Fermi envelope nxF, 
neglecting in particular the density hole due to the vortex 
line [3^ . Then one real izes that n^F^ A r, having a van- 
ishing divergence, see ( A42 ), can be written as the curl 
of some vectorial field B. One finds inside the Thomas- 
Fermi condensate ||l^ 



B 



Ng 2 



2mtj^ 



TF 



(A54) 



and one takes B ^ out of the Thomas-Fermi conden- 



sate. One then uses (A38) and (A36) to obtain 



cnrot 
phase 



TifiNgVl 



dsf(s) • e;n|p(ro(s)) (A55) 



where is the unit vector along z. Note that t (s) • = 
cosq;(s) where a is the angle between the vortex line 
and the rotation axis z, so that the integration element 
ds cos a{s) is simply dz, the length of the projection 
of the vortex line along z. From the estimate nxF ~ 
1/R\Rz ^ fJ-o/Ng and for a vortex length ~ R^ inside 
the Thomas- Fermi envelope one gets -Ephasc ~ —hQ, as 
expected from (L^) ~ h for a sin gle vortex configuration. 
This is in the energy range ( [A2| ) for a rotation frequency 
of the order of (|Al|). 



3. In terms of a 2D energy functional 

We now reexpress the 3D vortex line energy functional 
that we have obtained in terms of the energy functional 
of a 2D condensate with a vortex. Physically we asso- 
ciate a 2D fictitious condensate to each horizontal slice 
of the 3D condensate. Each 2D condensate is stored in 
a isotropic harmonic potential ma;^(a;2 -|- y'^)/2 with a 
rotation frequency and has the same number N of 
particles as the 3D condensate. The fictitious 2D con- 
densate at elevation z has a vortex at a position given 
by the intersection of the 3D vortex line with the hor- 
izontal plane of elevation z. The 2D condensate has a 
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Thomas-Fermi chemical potential 



stant of the 2D condensates: 



(A56) 



.92D (z) 



(ASS) 



where fio is the 3D vortex free Thomas-Fermi chemical 
potential. The Thomas-Fermi radius of the 2D conden- 
sate therefore coincides with the one R±{z) of the 3D We can then relate the 2D Thomas-Fermi density of the 



condensate at elevation z: 



R±{z) 



'2^J■2I){z) 



1/2 



(A57) 



Using the Thomas-Fermi value of the chemical potential 
for a 2D condensate we get the effective 2D coupling con- 



2D condensate to the Thomas-Fermi 3D density of the 
true condensate: 



/io - U{f) = gnTF{r) = g2i){z)n2B{r±; z). (A59) 



J 



Collecting all the energy terms of this section and after lengthy calculations we produce the result 



E-E' 



TF 



ds [E^^°iro^is);zo{s)) + cosia{s))E^,'^iro^is);zois))] 

9 



(A60) 



where a{s) is the angle between the vortex line and the axis z. The first 2D energy functional term contains the 
kinetic energy and the harmonic plus interaction potential energy. For a vortex core inside the Thomas-Fermi radius 
of the condensate, that is for r± = rj_/R±{z) < 1, it is given by 



E2^°{r±;z) 



where the value of the constant C is 



Ai2D(z) 



5 + {i-.l) 



C = Ao+Ai-l + -log2 + 



du 



C + lOE 



fJ'2I){z) 

huj I 



(1 - rl) 



F^{u) - 1 



(A61) 



(A62) 



Thomas-Fermi condensate and only E^]^^^^ has a non 



For > 1 the vortex line at elevation z is out of the 

s-Fermi condensate 
vanishing contribution to i?2ir°- 



^2B 



{rr,z) 



2M2d(2) 



1 + (1 - ri) log ^ 



1 



(A63) 



The second energy functional term in ( A60| ) is the rota- 
tional energy; it is given for < 1 by 



E','^{r^;z)^-hn{l 



~2 \2 

r I 



(A64) 



and it vanishes out of the condensate, that is for > 1. 
If one uses the hyperbolic tangent estimate of Q for F, 
see ( IaTbI) , one gets C ~ 0.0884. The 2D energy func- 



tional -E§5"° -t- £'20 then coincides exactly with the one 



of equation (64) in |6). 



4. Improving the 2D energy functional 

By a more detailed analysis than in |^ of the proper- 
ties of the 2D energy functional E^^^ we have identified 
a pathology that lead to problems in the full 3D energy 



functional minimization: the derivative of E^^^{r^]z) 
with respect to r± presents a logarithmic singularity on 
the border of the Thomas- Fermi condensate rj_ = i?j^(z), 
that is it diverges logarithmically to -l-oo in — 1^ 
and to —00 in fj^ = 1+. As a consequence E2-^'^{r_\_]z) 
has a local minimum inside the Thomas-Fermi conden- 
sate, at a distance of the order of the healing length 
C2D = fi/{mH2Diz))^^^ from the border, see Figure |2y. 
This local minimum of energy is an artifact of the ap- 
proximations used in the appendix: it is found by a nu- 
merical minimization of the full Gross-Pitaevskii energy 
functional with imaginary time propagation, that a vor- 
tex is going out of the condensate to infinity in the ab- 
sence of trap rotation. 

This artifact was not relevant in the 2D case, one just 
had to restrict to the Thomas-Fermi limit /i2D S> huj± 
with vortex cores inside the Thomas-Fermi condensate 
(excluding a thin layer of thickness ~ ^20 around the 
boundary). In the 3D case this artifact cannot easily be 
avoided: the bent vortex necessarily tries to cross the 
boundary of the condensate, where it encounters the log- 
arithmic singularity and can remain trapped. This prob- 
lem is also worse in 3D because the 3D vortex line can 
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(b) 




explore the extremities of the cigar shaped condensate, 
where the local 2D chemical potential fJ,2B{z) can be of 
the order of hLu±; the local minimum artifact in 
then becomes very pronounced, see Figure ^o|b, which 
seriously prevents the vortex line from bending and leav- 
ing the condensate. 



Heuristically we perform a modification to the energy 
functional i?^ff " eliminating the spurious trapping of the 
vortex line in the absence of rotation. In the region out of 
the Thomas-Fermi condensate, rj_ > 1, we realized that 
E^^^ is well approximated by a small 1 /f± expansion of 
the logarithmic term in (A63) so that we take 



2a*2d(2) 



1 

27^ 



2f4 



(A65) 



FIG. 20: For a 2D condensate in a non-rotating harmonic 
trap, energy of a vortex as function of the vortex distance 
r± to the trap center, for a fixed chemical potential (a) 
/i2D = 8huj± and (b) /12d = hu)±. So lid lin e: th e orig inal en- 
ergy functional of [p|, see equations (A61) and (A63), which 



presents a spurious local minimum inside the Thomas-Fermi 
condensate. Dashed line: the 77-regularized energy functional 
heuristically proposed here. 



In the inner Thomas-Fermi region, rj_<l, we eliminate 
the logarithmic singularity of the derivative at the inner 
border of the condensate by addin g a c onstant term t] to 
the argument of the logarithm of (A61): 
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El 



2D 



M2d(2) 



2 ri 
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1-r 



C + log 
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r 



These definitions lead to a continuous function E^]^^- 
The value of tj is adjusted to further ensure continuity of 
the derivatives of E^^'^ in f_L = 1: 



o3/4-C 



1.938. 



(A67) 



As seen on Figure ^ the regularized function has 
no local minimum, even for a chemical potential smaller 
than hcL!±. 

In 2D the physical predictions derived from i?^ff" 
slightly differ in the Thomas-Fermi regime from the ones 
from -E^ff". For example the critical rotation frequency 
il^^ such that the single vortex configuration has the 



same energy as the vortex free configuration is, from 
pa=o. 



^2D 



2D 



— ^log 

Ai2D 



=C+l/2 



M2D 
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^2D^Mlog 



(A69) 



whereas the original one in is 

^c+1/2 f mi 

M2D L" V^'^-L 

In the Thomas-Fermi limit f2^^ and fl^^ differ by a term 
scaling as rjti'^ u)^ / fi^D which is beyond accuracy of the en- 
ergy functional derivation of the present appendix. The 
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stabilization rotation frequency ^^^i, such that the single 
vortex core is a local minimum of energy is also modified: 



2D 
stab 



(A70) 



but this varies only to second order in rihuj±/ii2D- 



APPENDIX B: ENERGY OF THE VORTEX LINE 
SEGMENTS FAR FROM THE CONDENSATE 

We give here the behaviour of the field B defined by 
( [A44| ) far from the Thomas-Fermi condensate, that is at 
a distance r much larger than the Thomas-Fermi radii 
of the condensate. This problem is formally equiva- 
lent to the calculation in magnetostatics of the magnetic 
field B very far from a localized distribution of current 
j = nv/ iiQ where /iq is the magnetic permeability of vac- 
uum. This calculation is performed for example in 
From the property divj = one finds that the leading 
term in the l/r expansion corresponds to a dipolar field 



created by the magnetic dipole moment M of the current 
distribution: 



B{f) 



Mo 



3u{u- M] 



M 



(Bl) 



where u — r/r and the moment M is here proportional 
to the mean angular momentum per particle in the con- 
densate: 



1 



M = - / (Tf' r' A j (f ') = 



2miio 



(L). 



(B2) 



As a consequence B tends to zero as at infinity, so 
that the tot al flu x oi v A B vanishes at infinity, as it was 
assumed in ( A38 ). When the condensate wavefunction is 
symmetric under reflection with respect to the xy plane, 
as expected for a planar bent vortex contained in the xz 
plane, the x and y component of (L) vanish and we get 
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